5.1.1 栅格数据的基本处理

5.1.2.1 投影栅格

5.1.2.1.1 简述R-gis坐标系

## 定义

投影：
## EPGS:https://epsg.org/search/map
## PROJ4:http://spatialreference.org/
## 坐标系主要包含PROJ.4/EPSG两种形式，前者是定义详细的参数(R:raster/rgdal)，后者是直接查询对应定义的代码即可(R:sf)

## R :sp:CRS ## 定义坐标系：raster::CRS() 参数说明：
## 举例：proj4字符串投影定义：
+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## + proj = utm：投影为UTM，UTM有多个区域
## + zone = 11：区域是11，这是美国西海岸的区域
## +datum= WGS84：基准WGS84（基准是指投影中使用的坐标系的0,0基准）
## + units = m：坐标单位为METERS
## + ellps = WGS84：数据的椭球（如何计算地球的圆度）为WGS84
## 每个proj4字符串的最后一部分是+towgs84=0,0,0 。如果需要datum转换，则使用此转换因子(## 类似于七参数转换，这里是三参数转换)

## 分别使用proj4和ESPG定义地理投影：
## proj4:
a_crs_object <- crs("+proj=longlat +datum=WGS84 +no_defs")
class(a_crs_object)
## EPSG:
a_crs_object_epsg <- crs("+init=epsg:4326")
class(a_crs_object_epsg)

## 注意下载使用的数据一般为wgs84格式，但是谷歌地球使用的正常投影为Web墨卡托现在的正式官方代号 EPSG:3857；如果不能使用中国地球的投影面积，可以使用这个用于计算面积；
#

## 查看栅格的特征：
##查看单变量所有属性：
bio1 <- raster("data/bioclim/bio1.bil")
bio1 
nrow(bio1) 查看变量的行数；
extent(bio1)查看变量的范围;
crs(bio1) 查看变量的参考系
res(bio1) 查看变量的分辨率

5.1.2.1.2 栅格投影与投影转换

## 方法1：sp::spTRransform(x,crs())/sf::st_transform()
## 这个转换是对shp文件进行转换；
## sp::spTRransform(x,crs())
crs.laea <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs") 
locs.laea <- spTransform(locs, crs.laea)  
## sf::st_transform()
s.sf.gcs <- st_transform(s.sf, "+proj=longlat +datum=WGS84")
st_crs(s.sf.gcs)

## 方法2：raster::projection(x)/projectRaster(x,crs())
## raster::projection(x) <- value
## raster::projectRaster(x,crs())
crs.geo  ##包含CRS定义坐标系；
projection(tmin1.c) <- crs.geo
tmin1.proj <- projectRaster(tmin1.c, crs = crs.geo) 

## 方法3：raster::crs(x, ...) <- value
crs(occ_final) <- myCRS1 ## myCRS1是一个定义的CRS；


## 方法4：
library(sf)
library(tidyverse)
library(USAboundaries)

usa <- us_boundaries(type="state", resolution = "low") %>% 
  filter(!state_abbr %in% c("PR", "AK", "HI"))  # lower 48 only

wgs84 <- usa %>% 
  st_transform(4326) # WGS84 - good default

albers <- usa %>% 
  st_transform(5070) # Albers - a popular choice for the lower 48

ggplot() +
  geom_sf(data = wgs84, fill = NA)

ggplot() +
  geom_sf(data = albers, fill = NA)


## 方法5：常用地理坐标系转投影坐标系：
## 投影坐标系：
crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")  
proj4string(tiffs) <- crs.geo  

# 这里的工作就是将gs84转为下面的坐标系即可；
GCS_WGS_1984
WKID: 4326 权限: EPSG

## 参见wgs84转投影坐标系：
## wgs对应的投影是：
至今，EPSG:3857（WGS 84 / Pseudo-Mercator） 代号是web墨卡托的正式代号。
cord.UTM <- spTransform(cord.dec, CRS("+init=epsg:32647"))

EPSG4326：Web墨卡托投影后的平面地图，但仍然使用WGS84的经度、纬度表示坐标；
EPSG3857：Web墨卡托投影后的平面地图，坐标单位为米。

5.1.2.2 构筑范围/缓冲区/掩膜

##加载R包
library(raster)
##构筑范围（方法1：基于raster）
model.extent<-extent(min(rattler$lon)-10,max(rattler$lon)+10,min(rattler$lat)-10,max(rattler$lat)+10)
## 构筑范围(方法2:基于dismo)：
library(dismo)
> r <- raster(acg) ##注意这里的 acg是一个点文件shp；
> # set the resolution of the cells to (for example) 1 degree
> res(r) <- 1
> # expand (extend) the extent of the RasterLayer a little
> r <- extend(r, extent(r)+1)

## 根据分布点构建研究范围的若干方式
## 以下结果全部生成的shp文件
## 在使用时最好先调成对应的参考坐标系格式；
crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")  
proj4string(tiffs) <- crs.geo  

# 第一种方式：构建研究范围，一般以5度为限制；
## 注意这里定义了一个boundary的函数_frange;
frange <- function(x){
  occs <- x
  xmin <- min(occs$longitude)
  xmax <- max(occs$longitude)
  ymin <- min(occs$latitude)
  ymax <- max(occs$latitude)
  bb <- matrix(c(xmin-5, xmin-5, xmax+5, xmax+5, xmin-5, ymin-5, ymax+5, ymax+5, ymin, ymin), ncol=2)
  bgExt <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(bb)), 1)))
  return(bgExt)
}

# 第二种方式：构建最小面积buffer：MCP:
sp::coordinates(occs) <- ~ longitude + latitude
bg_mcp <- mcp(occs)
bg_mcp <- rgeos::gBuffer(bg_mcp, width = 0.5)


# 第三种方式：构建最小分布点周围buffer尺度；
occs <- xh_as[,2:3]
xmin <- min(occs$longitude)
xmax <- max(occs$longitude)
ymin <- min(occs$latitude)
ymax <- max(occs$latitude)
bb <- matrix(c(xmin, xmin, xmax, xmax, xmin, ymin, ymax, ymax, ymin, ymin), ncol=2)
bgExt <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(bb)), 1)))
# 增加buffer:0.5度~50km etc
bgExt <- rgeos::gBuffer(bgExt, width = 0.5)

# 第四种 划定研究范围：
# define circles with a radius of 50 km around the subsampled points
x = circles(subs[,c("longitude","latitude")], d=50000, lonlat=T)

## 使用raster包的crop掩膜：
modelEnv=crop(currentEnv,model.extent) ##currentEnv是总体环境图层；
modelFutureEnv=crop(futureEnv, model.extent)
##使用raster包的mask掩膜；与crop类似，mask也可以用于多个栅格同时掩膜；
bio1_mask <- mask(bio1, occ_buffer)
plot(bio1_mask)
plot(occ_buffer,add=T)
plot(occ_final,add=T,col="blue")

## ARCGIS中获取栅格变量的方法，也即raster to polygon；
a、System Toolboxes --> 3D Analyst Tools --> Conversion --> From Raster --> Raster Domain

b、System Toolboxes --> Conversion Tools --> From Raster --> Raster to Polygon

c、利用镶嵌数据集Footprint图层的方法来获取

# 获取研究范围的一种形式：
## 利用栅格图层相乘的方式获取研究范围：
## 这是因为空白值相乘会得到零值，从而在原有栅格范围中去除掉；
null.cells=ndvi*elev2*hfp2 #la multiplicación de todos los mapas propaga los valores nulos
plot(null.cells) #thi

5.1.2.3 栅格重采样

####栅格重采样####
##使用raster包的projectRaster来定义重采样
##连续性变量使用bilinear ；##类型变量使用ngb(nearest neighbor)
##这里的单位随着地形坐标系和投影坐标系之间差异；
new = projcectRaster(ss,crs =crs(pj),res =1000,method =bilinear)

## 查看栅格的分辨率：
library(maptools)
data("wrld_simpl")
plot(wrld_simpl)

###自定义重采样：
bio1 <- raster("data/bioclim/bio1.bil")
# define new resolution
newRaster <- raster( nrow= nrow(bio1)/10 , ncol= ncol(bio1)/10 )
# define the extent of the new coarser resolution raster
extent(newRaster) <- extent(bio1)
# fill the new layer with new values
newRaster <- resample(x=bio1,y=newRaster,method='bilinear')
# when viewing the new layer, we see that it appears coarser
plot(bio1)

5.1.2.4 栅格镶嵌

#Download two more tiles
srtm2 <- getData('SRTM', lon=13, lat=48)
srtm3 <- getData('SRTM', lon=9, lat=48)

#Mosaic/merge srtm tiles
srtmmosaic <- mosaic(srtm, srtm2, srtm3, fun=mean)

5.1.2.4 栅格重分类

##二值分类：
myLayer<- raster("data/bioclim/bio1.bil")
binaryMap <-   myLayer>= 100
plot(binaryMap)

##多值分类：
# values smaller than 0 becomes 0; 
# values smaller than 0 becomes 0; 
# values between 0 and 100 becomes 1; 
# values between 100 and 200 becomes 2;
# values larger than 200 becomes 3;

myMethod <- c(-Inf,  0, 0, # from, to, new value
              0,   100, 1,
              100, 200, 2,
              200, Inf, 3)
myLayer_classified <- reclassify(myLayer,rcl= myMethod)
plot(myLayer_classified)

## 利用栅格计算器重分类：
library(raster)
raster::calc
names <- list.files("C:/Users/admin/Desktop/tmean/",pattern = "tif",full.names = T)
for(i in 1:length(names)){
  assign(paste0("p",i),
         calc(raster(names[i]), function(x) { x[x<10] <- NA; return(x) }))

}

5.1.2.5 栅格计算

wet <- raster("data/bioclim/bio13.bil") # precipitation of wettest month
dry <- raster("data/bioclim/bio14.bil") # precipitation of driest month
# To calculate difference between these two rasters
diff <- wet - dry
names(diff) <- "diff"  ##注意这里的names命名很重要！！！
# To calculate the mean between the dry and wet rasters
twoLayers <- stack(wet,dry)
meanPPT1 <- calc(twoLayers,fun=mean)
names(meanPPT1) <- "mean"
# the following code gives the same results
meanPPT2 <-  (wet + dry)/2
names(meanPPT2) <- "mean"

5.1.2.5 栅格掩膜和裁剪

## 栅格裁剪raster:：crop()
## 注意栅格裁剪是范围裁剪，获取shp的边界极值后裁剪；
exten = CHN = getData("GDAM",CHN,level =0)
tif3 = crop(tif2 ,exten)

## 栅格掩膜raster::mask()
## 栅格掩膜，是实际shp面文件掩膜，获取shp的边界范围后掩膜
exten = CHN = getData("GDAM",CHN,level =0)
tif3 = mask(tif2 ,exten)

## 但是掩膜的环境分析过程，当掩膜体较大时效率交底：
## 解决方法是：
clim<-mask(crop(clim,bbox(bkg)),bkg)

5.1.2.6 栅格tif文件转asc文件

## 导出单个tif转asc
## 批量参见5.1.5.4 批量掩膜及asc导出；
library(raster)
f <- "path/to/downloaded/file.tif"
r <- raster(f)
writeRaster(r, "path/to/outfile.asc", format="ascii"，overwrite = T))

###  环境变量裁剪及掩膜(批量),导出tif和asc；
tifs <- list.files(path ="D:/XH/第三阶段_env/envsV3tif/",pattern="tif",full.names =T)
tiffs <- stack(tifs)

env_bgExt <- raster::mask(crop(tiffs$BIO17, bgExt),bgExt)
## 注意批量图层掩膜时，选择先用范围再裁剪的方法：
envs_bios <- mask(crop(tiffs,env_bgExt),env_bgExt)

d_list = as.list(envs_bios)
names <- names(envs_bios)
dir.create("D:/XH/fifth/sa_area")

export_masktif_asc <- function(x){
  names = names(x)
  output1 = paste0('D:/XH/fifth/sa_area/',names,'.tif')
  output2 = paste0('D:/XH/fifth/sa_area/',names,'.asc')
  writeRaster(x,output1,format ='GTiff' ,
              overwrite = T)
  writeRaster(x,output2,format ='ascii' ,
              overwrite = T)}

lapply(d_list,export_masktif_asc)

## 批量tif转asc：
## 列出文件名；
namesy <- list.files(path="./ssp126-202040",pattern = "tif")
namesy[1]
##多个环境因子转为asc后导出；
export_band<-function(x){
  names = namesy[x]
  output = paste0('C:/Users/chengshanmei/Desktop/ssp126-202040/',names)
  writeRaster(biocrops[[x]],output,format="ascii" ,
              overwrite = T)}
lapply(1:length(namesy),export_band)


## 并行方法：
## 
library(raster)
#### 环境变量构建tif及并行转asc：####
tifs <- list.files(pafenth ="D:/XH/莲子草/第三阶段_env/xhenvs2_5/",pattern="tif",full.names =T)
tiffs <- stack(tifs)
d_list = as.list(tiffs)
## 将环境变量转为asc格式：
dir.create('./env2.5_asc') 
getwd()
##需要注意下面的路径必须使用全拼路径：
library(parallel) #用于并行计算
library(snowfall)  # 载入snowfall包

names <- list.files(path ="D:/XH/莲子草/第三阶段_env/xhenvs2_5/",pattern="tif")
sfInit(parallel = TRUE, cpus = detectCores() - 4)
sfLibrary(raster) 
sfLibrary(base)
sfExport("names")
sfExport("d_list")

export_masktif_asc <- function(x){
  names = names(x)
  output2 = paste0('D:/Rtem/env2.5_asc/',names,'.asc')
  writeRaster(x,output2,format ='ascii' ,
              overwrite = T)}
sfLapply(d_list,export_masktif_asc)
sfStop()

5.1.2.6 栅格tif文件批量重命名

## 修改文件夹的内的文件名：
library(tidyverse)
newname<- list.files("E:/sjdata/F2_ENVS/aseutif/",pattern = "tif") %>% strsplit(.,".tif")
nw <- c()
for(i in 1:42){
  nw <- c(nw,newname[[i]][1])
}
nw

news <- paste0("E:/sjdata/F2_ENVS/aseutif/",nw,".tif")
file.rename(from = list.files(path="E:/sjdata/F2_ENVS/aseutif/",pattern = "tif",full.names = TRUE),
            to = news)

5.1.2.7 栅格tif文件主成分分析

## pca——tiff：pca提取栅格：
library(sp)
library(raster)
library(maps)
library(mapdata)
## install.packages("RStoolbox")
library(RStoolbox)

## 环境数据集：
tifpca <- list.files(path ="D:/XH/莲子草/第三阶段_env//",pattern="tif",full.names =T)
tiffpcas <- stack(tifpca)

pcamap<-rasterPCA(tiffpcas,spca=TRUE)
summary(pcamap$model) ## pc1+pc2的解释度只有77%；

ff <- unlist(pcamap$map)
writeRaster(ff$PC2,"./pc2.tif")
## 已经导出来了；基于这两个pca去做相关性分析；

# 方法2：进一步结合pearson相关分析结果，选择|r|<-.8以下的变量；

## 方法3：不考虑共线性转移的问题，直接使用专家意见，构建限制物种分布的关键数据集；